Solid-Liquid Interface Free Energy through Metadynamics simulations 
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The solid-liquid interface free energy 7^ is a key parameter controlling nucleation and growth 
during solidification and other phenomena. There are intrinsic difficulties in obtaining accurate 
experimental values, and the previous approaches to compute 7 S ; with atomistic simulations are 
computationally demanding. We propose a new approach, which is to obtain y„i from a free energy 
map of the phase transition reconstructed by metadynamics. We apply this to the benchmark case 
of a Lennard- Jones potential and the results confirm the most reliable data obtained previously. We 
demonstrate several advantages of our new approach: it is simple to implement, robust and free of 
hysteresis problems, it allows a rigorous and unbiased estimate of the statistical uncertainty and it 
returns a good estimate of of the thermodynamic limit with system sizes of a just a few hundred 
atoms. It is therefore attractive for using with more realistic and specific models of interatomic 
forces. 

PACS numbers: 64.10.+h,31.15.xv 



I. INTRODUCTION 

Many important phenomena occurring in first order 
phase transformations, such as nucleation and growth, 
are controlled by interfacial properties. In the theory 
of solidification, one such property is the solid-liquid in- 
terfacial free energy 7^. This parameter controls the 
barrier for nucleation of a solid in an undercooled liq- 
uid and the transitions between planar, cellular and den- 
dritic growth regimes in metals, which in turn governs 
their final microstructureM Despite its importance for 
both theoretical models and practical applications, accu- 
rate data for the value of 7^ are not known even for the 
case of simple elements. There are indeed few experi- 
mental techniques aimed at measuring this quantity (for 
a comprehensive review see RefPj) and their application 
is complicated by the very strict control on all experimen- 
tal parameters that must be achieved to obtain accurate 
data. One such method for example involves recover- 
ing 'ysi indirectly from nucleation-rate measurements^. 
In this case, large uncertainties in the measured values 
arise from the possible occurrence of heterogeneous nu- 
cleation from very low-concentration impurities. Reliable 
theoretical values would therefore be very useful. 

Several methods have been developed to calculate 
jsi from in-silico experiments with molecular dynamics, 
where complete control of the "experimental" variables 
is achievable. These methods are the Capillary Fluctua- 
tion method (CFMj^j, different sorts of so-called "cleav- 
ing" methods (CM^ an d a Classical Nucleation The- 
ory (CNT) approaclP. In CFMs the fluctuation spec- 
trum of the interface height is related to the interfacial 
stiffness 7 S ; (9) + 7"; (9) (where the second derivative is 
taken with respect to an angle 9 defining the crystallo- 
graphic orientation of the surface) through which 7^; can 
be recovered by calculating 7^; + 7^ for different inter- 



face orientations and fitting the results to an expansion 
of 7 S ; in kubic harmonics^. In CMs, as the name sug- 
gests, bulk solid and liquid phases are separately cleaved 
and the different phases are joined to form an interface. 
In this way, 7 S ; is recovered by measuring the work done 
during the process. Finally, in the CNT approach, crys- 
talline nuclei of different sizes are inserted into a super- 
cooled liquid and some orientational average of 7 S ; is re- 
covered by measuring the radius of the critical nucleus 
R* and inserting its value in the classical nucleation the- 
ory equation relating R* and 7^ (see for example RefP, 
page 46) . We refer the interested reader to the literature 
for details of these calculations. Successful applications 
of the aforementioned methods hav e bee n reported for 
model systems such a s hard sphere d 4 I 5 I 9 I and Lennard- 
Jones potential d 6 * 10 * 11 ! as well a s mo re realistic semiem- 
pirical and quantum-mechanicaPE^^ based Embedded 
AtoirP^ and Stillinger-Weber^ potentials. 

The CFM and CNT are derived with macroscale ap- 
proximations and thus require large simulation supercells 
of about 10 5 atoms to be applicable and give accurate 
results. Cleavage methods require somewhat smaller su- 
percells (« 10 3 — 10 4 atoms) but are prone to the error in- 
troduced if the sequence of simulations is not completely 
reversible. A dramatic example would be the complete 
solidification of the liquid while joining it to the solid 
due to a large relative fluctuation in the position of the 
interface^. The simulation supercell must contain a rel- 
atively large area of interface in order to avoid the oc- 
currence of these events. Moreover, to compute j s i ac- 
curately and efficiently, one has to use a cleaving po- 
tential which mimics accurately the interactions between 
the system's particles^. This must be built in an ad hoc 
way for every system and can become cumbersome when 
complex many-body interactions have to be taken into 
account such as for example in ab-initio-based calcula- 
tions. 
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These shortcomings become particularly troublesome 
if one consider that interface free energies are very sensi- 
tive to the details of the empirical potential; for instance, 
different paramctcrizations of EAM potentials yield val- 
ues of j s i which vary by as much as 30% 19 . In order to 
capture the complex bonding and the unusual local en- 
vironments present at the solid-liquid interface, and to 
capture accurately the anisotropy of crystalline surface 
energies, one must consider more sophisticated models, 
which reproduce more closely the first-principles total en- 
ergy. 

In the present paper, we discuss a novel technique 
to compute 7^ which aims at being robust, efficient 
and transferable, and which is a promising candidate 
to extend the scope of interfacial energy calculations 
to more complex potentials than previously treated. 
Briefly, our method reconstructs a coars e-grai ned free 
energy surface (FES) using metadynamics^EU. Such a 
FES maps out the transition from a single phase to the 
space of configurations where two phases coexist. The 
minimum difference in Gibbs free energy between these 
two regions at the solid-liquid equilibrium temperature 
is an excess free energy G xs , which is equivalent to the 
interface free energy j s i multiplied by the area A of the 
interface. 

The remainder of this paper is organized as follows. In 
Section II we present the thermodynamic basis and the 
details of the method. In Section III we describe the com- 
putational details of our simulations. In Section IV we 
show our results for a simple Lennard- Jones system and 
critically discuss them in comparison with other available 
methods. We also speculate on the possibility of imple- 
menting this approach together with ab-initio molecular 
dynamics. Finally, we summarize our main results. 



at T m is identical, one can write the overall Gibbs free 
energy as 

G s \i(P,T m ) = fi(P,T m )N + G xs , (2) 

where an excess energy term associated with the interface 
has been introduced. Such a term will be extensive with 
respect to the area of the interface, and we can write it 
as the product of the surface area A and an interface free 
energy j sh i.e. G xs = Aj sl . 

The most direct approach to the computation of j s i 
is clearly to calculate the free-energy difference between 
the bulk phases and the configurations in which planar 
interfaces are present, as described by Eqs. [T] and [2] re- 
spectively. We will obtain this free-energy difference by 
means of metadynamics simulations, as described in the 
next section. 



B. Free energy differences from metadynamics 

The use of metadynamics for reconstructing free- 
energy landscapes has been the subject of many papers 
and we refer the reader to the excellent review by Laio 
and Gervasio and references thereirP^l, while we only 
briefly sketch the main ideas here. Metadynamics is a 
simulation technique based on non-equilibrium molecu- 
lar dynamics, which is designed to reconstruct a coarse- 
grained free energy surface (FES) in the space of one 
or more collective variables {si} that describe the state 
of the system. Metadynamics reconstructs the FES by 
adding a bias potential in the form of a Gaussian cen- 
tered at a specific point in the Collective Variable (CV) 
space each time that point is visited. The mathematical 
form of the bias potentials is given by 



II. METHODOLOGICAL DETAILS 

A. Thermodynamic basis 

We consider a homogeneous solid or liquid system of N 
atoms, located in a periodically repeated supercell within 
an infinite system, at a pressure P and temperature T. 
Its Gibbs free energy G can be written as 



G s{l) (P,T) = » s{l) (P,T)N 



(1) 



where fi s rn is the chemical potential of atoms in the 
solid (liquid) phase. At the melting temperature T m , 
the chemical potentials in the two phases are equal 

fi s (P, T m ) = m (P, T m ) = M (P, T m ) . 

There exists a second state of the same system at the 
melting temperature, in which solid and liquid phases 
coexist, separated by macroscopically planar interfaces 
that are naturally fluctuating on the atomic scale. Since 
the chemical potential in the solid and liquid bulk phases 
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which in the discrete version needed to implement the 
algorithm for computations becomes 
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Here r is the inverse of the frequency of deposition of 
the Gaussians, and N = i/r is the number of Gaussians 
accumulated up to time t in the simulation, w is the 
deposition rate of the Gaussian functions and a their 
width. 

The bias discourages the trajectory from remaining in- 
definitely in the same region of the CV space, effectively 
pushing the system towards the lowest-lying free-energy 
barrier. Once all the relevant free energy minima have 
been levelled by the bias (see Figure I]), the system be- 
comes completely diffusive and wanders freely through 
all the possible states. 
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Figure 1: (color online) Schematic representation of the flat- 
tening of the effective FES by the repulsive bias of Eq. Q, 
during a met adynamics simulation in a one-dimensional col- 
lective variable space. We show the underlying FES G(s) and 
the bias accumulated at chosen times (arbitrary units). For 
a sufficiently long simulation, G(s) + V(s, t) — > constant, so 
that one can obtain an accurate estimate of the free-energy 
surface simply by taking the negative of the bias. 



At this stage of the simulation the accumulated bias is 
equal to the negative of the free energy of the real sys- 
tem plus an additive constant (for a detailed analysis see 
RefJ22l). However, there are two limitations in this orig- 
inal form of metadynamics. First of all it is not clear 
when a metadynamics simulation should be stopped, i.e. 
when the bias has effectively compensated the underly- 
ing free energy. Moreover, even at this point, the effec- 
tive FES will have a residual roughness of the order of 
the height of each individual Gaussian (wt in equation 
|4|. In order to resolve these issues, the so-called "well- 
tempered" metadynamics methocpf has been proposed 
recently by Barducci et al., and we use this specific type 
of metadynamics in our simulations. The idea behind 
well-tempered metadynamics is to gradually reduce the 
height of the deposited Gaussians, at a rate determined 
by the magnitude of the bias already present. The ex- 
pression analogous to ^ reads 
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The parameter AT controls how quickly the deposition 
rate is reduced. Once the simulation approaches conver- 
gence, the collective variables space will be sampled with 
a probability distribution corresponding to an artificial 
temperature T + AT^I Hence, the final bias accumu- 
lated during a single simulation converges to 
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The true free energy can be recovered inverting equa- 
tion 

As in any free-energy calculation based on the map- 
ping of the configurations of the system to a coarse- 
grained collective-variable space, the definition of CVs 




Figure 2: (color online) Angular function c a (x) as defined in 
equation [8] The function is shown as a polar plot, centered 
on an fee lattice. c a has well-defined peaks in the directions 
of the nearest neighbors. 



that can effectively distinguish between relevant states, 
and describe reliably the natural transformation path is 
the first, and most important step. The primary require- 
ment is to distinguish the solid phase from the liquid. 
With this aim, we define for every atom an order pa- 
rameter </>, which depends on the relative position of its 
neighbors. The definition of <f> 

E«i °r (l x j ~ x i|) c a ( x j — x i) 

<M X ») = — r~nz — zrvi ■ ( 7 ) 
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contains an angular term c a to distinguish the different 
environments, and some radial cutoff functions c r which 
are useful to guarantee that is a continuous function 
of all its arguments. Note that the weighted sum of c a 
is normalized over the total coordination, so that <fi is 
relatively insensitive to fluctuations of the density. 

We define the angular function c a as a combination of 
polynomials in Cartesian coordinates, symmetry adapted 
to the cubic point group: 

c a (x) = [x A y 4 (l - z 4 / |x| 4 ) + yV (l - .t 4 / |x| 4 ) 

1 (8) 
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We have chosen Eq. ([8]) rather than more traditional 
parameters such as the so-called Q§ (see e.g r^ 2 ° j) ! for 
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Figure 3: (color online) Cutoff function used to define the 
regions A and B for the calculation of the two collective vari- 
ables. The function varies smoothly from to 1 so as to avoid 
discontinuities when atoms transit between the two regions. 



We finally define our CVs sa and sb by averaging the 
order parameters of the atoms located within region A 
and B, respectively: 
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where 
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This scaling has been chosen so that <j> = in a homoge- 
neous liquid and <f> — 1 in a perfect fee solid. 



III. COMPUTATIONAL DETAILS 



a number of reasons: c a has well-defined peaks for an 
fee environment (see Figure [2| , it is not rotationally in- 
variant (and will therefore enforce an orientation of the 
crystal consistent with the periodic boundaries) and it is 
relatively cheap to compute. It would possible to con- 
struct a different form of c a if one wanted to deal with 
a different crystal structure, and one simply has to ro- 
tate the function in order to specify a different crystal- 
lographic orientation of the surface. The application of 
a specialised, orientation-dependent order parameter is a 
key ingredient of our approach. 
The radial cutoff is defined as 
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where y = (r — ri)/(fo — ri). The polynomial part in 
Eq. (J9j) is simply a third order polynomial satisfying the 
constraints of continuity of cy(r) and its first derivative 
at r\ and tq. 

In order to study the formation of a solid-liquid inter- 
face, one must then distinguish configurations where the 
supercell is completely solid, completely liquid, or par- 
tially solid and partially liquid: in the latter case, at 
least two parallel interfaces will be present. For this pur- 
pose we divide the supercell, centered at the origin, into 
two regions: we assign to region A those atoms having 
\z\ < z, and to region B all the others (see Figure [3]). 
Note that we take z to be about one fourth of the super- 
cell length along z, and we keep it constant irrespective 
of the fluctuations of the supercell's size. This choice is 
not troublesome as long as the averages are properly nor- 
malized, so that the value of the CVs is independent of 
the number of atoms contained in each of the regions. 

Again, in order to obtain smoothly-varying CVs, we 
introduce a weight function. We use the same functional 
form introduced for the radial cutoff; namely, c z (x) = 
(Vd^l), setting n = z, r — z + Az in Eq. 



In order to evaluate our method, we decided to per- 
form the metadynamics calculations with a truncated 
Lennard- Jones potential, in the form used by Broughton 
and Gilmep22l. Such a simple potential is inexpensive 
and thoroughly studied, yet it can capture many of the 
relevant physical phenomena occurring in real systems. 
Its solid-liquid transition, an important ingredient of 
our approach, has been recently studied by Mastny and 
de PabkP3 through a modified Wang-Landau sampling 
technique^. Therefore, this system is ideal for a sys- 
tematic study with our method, comparing it with other 
techniques and benchmarking its performance as a func- 
tion of the parameters of the simulation. 

We will use Lennard-Jones units throughout the pa- 
per. The zero pressure coexistence temperature for this 
system has bee n rece ntly recalculated and we take it to 
be T m = 0.6188 111221 . Details of the phase dia gram can 
be found in refP^. We performed our simulations with a 
range of supercell sizes from 4x4x6 fee unit cells (384 
atoms) to 9x9x20 (6480 atoms). The supercells were ori- 
ented with fee [001] cell vectors parallel to the axes, with 
the longest side parallel to z, and were rescaled to a vol- 
ume consistent with the equilibrium density of the solid 
at the coexistence temperature. Atomic positions were 
then equilibrated at T m by performing a short molecu- 
lar dynamics simulation in the canonical (NVT) ensem- 
ble. This procedure was adopted in order to generate the 
starting configurations for the subsequent metadynamics 
simulations, which we perform instead in the isothermal- 
isobaric (NPT) ensemble. The timestep for the integra- 
tion of the equations of motion, performed with a velocity 
Verlet algorithm^, was 0.004. This choice gave negligible 
drift of the conserved quantity in all our simulations. 

In order to perform constant pressure simulations, 
variable-cell dynamics is implemented using a Langevin- 
piston barostatP^ and a friction of jb = 2 ps _1 . 

The presence of an interface calls for particular atten- 
tion when performing constant-pressure simulations. In 
particular, one has to deal with the change in density 
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that occurs when a portion of the supercell melts, at the 
same time considering fluctuations in the a;y-plane. If the 
xy-plane parameters of the supercell are fixed, the fluctu- 
ations in the solid will be frustrated; on the other hand, 
if those parameters are left free to vary, one will witness a 
spurious shrinking of the dimensions in the xy-plane due 
to surface tension whenever an interface is present. In 
both cases one can in principle observe a strain-related 
free energy contribution. However, this problem would 
disappear in the thermodynamic limit, hence one can just 
make the choice that is more computationally convenient, 
and consider the resulting error as another finite-size ef- 
fect, which can be controlled by comparing the results 
with different supercell sizes. We decided to let only the 
z component free to fluctuate. In this way, the change 
of volume occurring as the fraction of liquid and solid 
phases changes can be accommodated, and even in case 
of complete melting the xy dimensions remain commen- 
surate with a strain-free solid, which makes it simpler for 
the system to freeze again into a defect-free solid. 

Temperature control is extremely important in meta- 
dynamics simulations, since the increase of the biasing 
potential creates a continuous supply of energy to the sys- 
tem, which must nevertheless be held close to equilibrium 
in order to sample the free energy correctly. A strong lo- 
cal thermostat is needed, but at the same time one must 
avoid overdamping, which drastically reduces the diffu- 
sion coefficient and hence the sampling of slow, collective 
modes. We therefore use a colored-noise thermostatP^U 
fitted to provide efficient sampling over a broad range 
of frequencies, corresponding to vibrational periods be- 
tween 0.1 and 10 3 Lennard- Jones time units. 

The metadynamics we used for the different simula- 
tions are reported in Table [l] This Table also includes 
data for a number of tests using a single CV, which we 
describe later. We performed tests with other choices of 
these parameters spanning about an order of magnitude 
and no statistically significant changes were observed 
in the calculated value of 7^. The values reported, 
however, resulted in the best statistical uncertainty 
in the final free-energy estimate. The simulations were 
performed using the DL_POLY code (version 2.18p^), 
patched to perform metadynamics using the PLUMED 39 
cross-platform plugin, which greatly reduces the imple- 
mentation burden by providing a convenient framework 
for introducing new collective variables. 

When performing simulations at T > OK, thermal fluc- 
tuations induce some disorder in the solid and the scaled 
order parameter <p deviates from the value predicted for 
the ideal fee crystal. In figure [4] we report the time- 
averaged order parameter (d>j and its fluctuations for a 
single atom in the bulk phases. At the coexistence tem- 
perature, the average for the liquid oscillates around zero, 
while for the solid {£) ps 0.83. We note that even for an 
individual atom the order parameter can distinguish very 
clearly between the two phases at the melting tempera- 
ture. 



# atoms (cell) r 


1 + ^ 


WT 


SI (2D) 2352 (7 x 7 x 12) 4 


60 


0.115 


S2 (ID) 384 (4 x 4 x 6) 4 


40 


0.037 


S3 (ID) 512 (4x4x8) 4 


40 


0.037 


S4 (ID) 768 (4 x 4 x 12) 4 


40 


0.037 


S5 (ID) 1024 (4 x 4 x 16) 4 


40 


0.037 


S6 (ID) 1280 (4 x 4 x 20) 4 


40 


0.037 


S7 (ID) 1200 (5 x 5 x 12) 4 


60 


0.058 


S8 (ID) 2352 (7 x 7 x 12) 4 


120 


0.115 


S9 (ID) 6480 (9 x 9 x 20) 4 


205 


0.191 



Table I: Metadynamics parameters for different simulations 
(1 and 2 dimensional) and different supercell sizes. AT has 
been chosen such that k(AT + T) « -y 3 iA for every size. An 
order-of-magnitude estimate of 7^ suffices for this purpose, r 
was chosen so as to observe the first solid-liquid transition at 
about half of the total simulation time, and a was set to 0.03, 
which is of the order of the thermal fluctuations of the CVs 
in an unbiased simulation. 
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Figure 4: (color online) Time-averaged value of the order pa- 
rameter <j> for an individual atom in the bulk solid and liquid 
phases, as evaluated for the LJ system at different temper- 
atures and across the solid-liquid transition. The bounding 
lines delimit the range one standard deviation above and be- 
low the mean value. Dashed lines correspond to the values 
of the order parameter for metastable solid and liquid. Note 
that the parameters r\ and ro for the radial cutoff (9| are 
scaled from the values used at T m according to the changes 
of the equilibrium density. 



The parameters entering the radial cutoff function c r 
have been chosen to be r = 1.5 and r\ = 1.2, so as to 
encompass the typical first-neighbour distances in both 
Lennard- Jones solid and liquid at T — T m . In order to 
prevent sa and sb from visiting irrelevant configurations, 
corresponding to an order parameter far from its mean 
value in either liquid or solid, we have applied a lower 
and upper wall on both CVs^l in the form 



V wa u(s) = k 



S S limit 



with k = 50, e = 0.01, n = 4 and sn 



(12) 



-0.15 and 
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0.95 for the lower and upper wall respectively, which in- 
troduces a restraining potential for sa,sb < —0.15 and 
sajSb > 0.9. At the same time, V wa u is set to be zero 
inside this interval and thus cannot interfere with the 
dynamics of the system in this region of CV space. 



IV. RESULTS AND DISCUSSION 
A. Qualitative analysis of the FES 

Ideally the FES should be symmetric about the line 
sa = sb- Moreover, as calculations are performed at 
T m , one should observe the occurrence of two minima 
with the same free energy, at the values of the CVs cor- 
responding to the single bulk phases (either solid or liq- 
uid). The expected behavior is clearly exhibited by the 
calculated FES, reported in Figure [5] for a 7x7x12 super- 
cell, where we show the free-energy landscape together 
with some snapshots corresponding to different values of 
the CVs. The combination of CVs s = {sa + sb)/2 cor- 
responds roughly to the average of <f> over the whole box, 
and distinguishes between configurations with different 
proportions of solid and liquid phases. It can be used 
as a convenient reaction coordinate (see Figure [5jf ) for 
the FES projected on s). As expected, two wells occur 
with minima at the complete solid and liquid states, sep- 
arated by a rather flat region, whose height above the 
two minima corresponds to the interfacial free energy. 

The fact that the two wells should have the same depth 
can be used to check that the simulation temperature is 
indeed the melting temperature. This is a significant ad- 
vantage of our method, since knowledge of the melting 
temperature is a prerequisite of all the different meth- 
ods described in the literature to calculate 7 s j . Both the 
CFM and CM, being based on coexistence simulations, 
could suffer from major irreversible changes leading to 
complete solidification/melting of the simulation cell if 
not performed at the correct temperature, and the data 
gathered during the simulation would not serve its pur- 
pose. Also the CNT method needs the correct value for 
the melting temperature in order to calculate the nucle- 
ation barrier from which j s i is recovered. Our method, 
by contrast, is still effective even if the simulation tem- 
perature is slightly off the actual T m . Clearly, an error 
will be introduced, since "/ s i is estimated from equation 
[2] which is satisfied exactly only at T = T m . However, 
our method is very robust, in the sense that it tells us 
both whether such an error occurs and give us the sign 
and an estimate of the magnitude of the correction. We 
will discuss this issue further when addressing finite-size 
effects. 

The two-dimensional FES is rather constant along 
the line of points equidistant between the two wells, 
in the direction of the orthogonal combination of CVs 
s = (sa — Ss)/2, since this variable describes the posi- 
tion of the two phases with respect to the partitioning of 
the cell along z (see snapshots a, b and e in Figure [5]) . 



As we will comment further below, our CVs are effec- 
tive because there exist no metastable phases of the LJ 
potential correspond to values of the order parameter </> 
between 4>i and <f> s , and this is likely to be the case for any 
potentials describing a fee crystal. Hence, when moving 
away from the perfect liquid(solid) bulk value, any ho- 
mogeneous variation of the order parameter would be 
too energetically costly, and the system instead induces 
some order (disorder) in the form of small clusters, slowly 
increasing the free energy (region in orange in Figure [5| . 
Because of the elongated aspect ratio of the supercell, as 
soon as enough liquid(solid) phase is present, the most 
favourable configuration corresponds to the presence of 
two interfaces perpendicular to the z axis. The fact that 
we observe a solid<-> liquid transition via the growth of an 
individual cluster suggests that a very similar approach 
can be used to study the nucleation process itself. This 
idea will be explored in future work. 

As the time needed by metadynamics to reconstruct 
the FES is an exponential function of the dimensionality 
of the coarse grained space, one might wonder if it is pos- 
sible to speed up calculations by using a single CV. Wc 
explored this possibility as follows. Rather than using 
s = (sa + ss)/2, we kept a two-dimensional description, 
but performed metadynamics on sa alone, while atoms in 
region B are constrained in order to maintain this region 
of the supercell in a solid state (i.e we apply a restrain- 
ing lower wall potential which is a function of sb, and 
introduces a penalty in the enthalpy whenever sb devi- 
ates too much from <f> s . The values of the parameters 
entering in the restraining potential, whose form is de- 
scribed in Section |Hl} are k = 50, e = 0.01, n = 4 and 
s Hmit — 0.7. This forces region B to remain solid, while 
region A can sample both solid and liquid phases. In 
this case, the FES should show a minimum for sa ~ 4> s 
where the supercell is completely in a solid phase and a 
maximum where sa ~ 4>i i.e. when two interfaces are 
present. Again the difference between these two values 
is the interface excess energy. Figure [8] (inset (a)) shows 
the ID FES reconstructed in this way at different simu- 
lation times. The use of a single CV does not have any 
adverse influence on the calculated value of 7 S ;, as we will 



show in Section IV B and the use of this simpler form of 



metadynamics is fully justified for our purposes. 

A necessary condition for metadynamics to reconstruct 
the coarse-grained free energy of the system in a meaning- 
ful way is that all the important states and the barriers 
between them are effectively reached many times during 
the simulation. Moreover, one has to make sure that 
quasi-equilibrium conditions hold, which can be moni- 
tored by checking temperature and structural relaxation 
of the system. 

In order to check that the system effectively performs 
many transitions between the single-phase and the two- 
phase states, we verified that the CVs oscillate several 
time between their value in the liquid and solid phases. 
Moreover, we also visually check that the system actu- 
ally performs these transitions by printing snapshots of 
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Figure 5: (color online) 2D FES reconstructed by well-tempered metadynamics, together with selected snapshots of configura- 
tions obtained during the simulation. Atoms participating in sa are colored in orange, those in sg in blue, and the region of 
CV space corresponding to each snapshot is marked. The negative peaks in the FES clearly correspond to the two single-phase 
regions. They are separated by a very wide plateau, corresponding to the presence of well-defined interfaces between solid 
and liquid phases at various different positions relative to the A/B partition (insets (a), (b), (e)). In inset (f) we report the 
projection of the FES along the single CV s — (sa + ss)/2. 
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Figure 6: (color online) Radial pair correlation function g(r) 
for the liquid in the presence of an interface during our meta- 
dynamics simulations (red, dashed) and for a normal molec- 
ular dynamics simulation of a bulk liquid (blue). It is clear 
that the two curves are very similar, thus ruling out quan- 
titatively the presence of metastable structures during our 
simulation. In the inset, the kinetic energy distribution dur- 
ing our simulation is plotted in comparison to that expected 
for a canonical ensemble at T = T m . Again the two are very 
close demonstrating that the metadynamic bias does not in- 
duce any systematic deviation from the correct ensemble, and 
that quasi-equilibrium conditions hold. 



the atomic positions and visualising them using the VMD 
software^. Quasi-equilibrium conditions hold very accu- 
rately, as demonstrated from the inset (b) of Figure [6] 
where we show the velocity distribution function com- 
pared to its analytical equilibrium value. The radial pair 
correlation function g(r) of the liquid portion of inter- 
facial configurations (Figure [6| agrees well with the one 
computed for the bulk liquid in an unbiased run, which is 
a further confirmation that our simulation strategy does 
not introduce spurious structural effects. The g(r) dis- 
tribution is a sensitive measure of the short-range order 
present in the liquid, and any extra structuring would 
have been clearly detected as a shift of the peak positions 
or shapes, which does not happen here. The absence of 
artefacts has also been checked by visual inspection of 
snapshots of the atomic configurations along the meta- 
dynamics trajectory. 

A peculiar feature of our approach is that, at vari- 
ance with cleavage methods, the solid-liquid interface is 
created and "annihilated" several times during each sim- 
ulation, so that hysteresis should be much less of a con- 
cern. When the well-tempered bias is nearly converged, 
the systems diffuse on a flattened FES, and the morphol- 
ogy of the interface corresponds to the most favourable 
one from a free-energy point of view. As seen from Fig- 
ure [7J such a morphology includes a significant amount 



Figure 7: (colour online) A snapshot of the solid-liquid inter- 
face taken from the final part of a well-tempered metadynam- 
ics simulation. The scaled order parameter (p has been used 
to colour atoms. The atoms with a liquid-like configuration, 
with (j> < 0.45 have been hidden, the atoms with 4> > 0.65 
have been coloured in blue. Finally, atoms in intermediate 
configurations, with 0.45 < (f> < 0.65 have been made translu- 
cent. It is clear that - whatever threshold is used to ascertain 
the solid from the liquid state - the interface is not flat on the 
atomic scale. 

of roughness at the atomic scale. This might be expected 
from the observation that for the system under consider- 
ation the melting temperature is higher than the rough- 
ening transition temperature for the (100) surface. 

B. Analysis of accuracy and system-size effects 

Several terms contribute to the error in calculating a 
complex thermodynamic property such as 7 S ;. In actual 
applications of this method to a real substance one will 
be concerned with the accuracy of the total energy and 
force model, but this is not an issue in our present proof- 
of-principle case. However, there are still two major 
sources of error we must be concerned with here; namely, 
a statistical error stemming from insufficient ergodicity 
of the sampling (a finite sampling-time error) and the in- 
accuracies caused from insufficient size of the supercell. 
These finite-size errors introduce a lower bound on the 
acoustic vibrational frequencies, and most importantly 
might affect the structure of the liquid phase, introduc- 
ing spurious correlation that change the liquid entropy 
thus changing the melting temperature of the system (al- 
though we have seen no evidence for this in the pair cor- 
relation function reported above). Moreover, they could 
in principle induce a strain field in the solid and introduce 
interactions between the two interfaces. 

The finite-sampling error is readily gauged, by per- 
forming several independent runs and by checking how 
quickly the discrepancy between the reconstructed free- 
energies converges to zero. It is shown in Re fs 12112H th a t 
for simple models the error in the FES, after a short tran- 
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Figure 8: (colour online) Convergence of the FES (inset (a)) 
and its average error (inset (b)) with respect to time for our 
ID simulations. Ten simulations have been performed on a 
7 x 7 x 12 supercell, with the single-CV setup described in 
the text. The FESs in inset (a) are constructed by averaging 
equal-times biases of the independent runs, and the error- 
bars correspond to the standard deviation. In inset (b) we 
plot such an error, averaged between <f>i and (j> s , as a function 
of simulation time. The error is plotted on a log-log scale and 
the least-square linear fit shows that the angular coefficient is 
close to the theoretical value of —1/2 predicted for a simple 
Langevin model. 

sient phase when the free-energy basins are being filled, is 
expected to decay as the inverse square root of simulation 
time. 

It is reassuring to verify that this behaviour is also 
found in our system, as shown in Figure [8] This be- 
haviour is to be expected, because for long simula- 
tion times well-tempered metadynamics corresponds to a 
histogram-reweighting with a nearly perfect biasing po- 
tential. It means also that rather than running a very 
long simulation, one can with equal machine efficiency 
perform several, shorter, independent runs, with great 
advantages from the point of view of parallelisation. 

In order to calculate the value of 7 S ; from the recon- 
structed FES, we need to monitor in time the estimate 
of the excess free-energy due to the interface, 

G xs (t) = G s{l (t) - G s{l) (t) -> lsl A , (13) 

where we label as G s \i the estimate of the free energy for 
a configuration with a solid-liquid interface. As is rou- 
tinely done in conventional metadynamics simulations, 
we take as our best estimate of G xs the incremental av- 
erage over the final part of the trajectory, well after the 
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initial transient: 

G xs « — !— [ f G xs (t)dt. (14) 

We perform ten independent runs, and we can therefore 
compute an unbiased estimate of the overall statistical 
error. 

As we previously discussed, in the case of our 2D meta- 
dynamics, there are many points on the FES correspond- 
ing to coexisting solid and liquid phases, which have the 
same free energy (see Figure [5]); analogously, an extended 
plateau region is found in the ID setup. Therefore, any 
point in these regions would be a valid choice for evaluat- 
ing the interfacial free energy, provided that these regions 
are indeed flat. This brings us to the discussion of finite- 
size errors. In fact, at least for a simple, short-ranged 
potential such as Lennard- Jones, the greatest concern is 
the interaction of the two interfaces along z, mediated 
by the elastic strain field in the solid portion and by the 
altered structure of the liquid in close proximity to the 
solid/liquid boundary. Such effects are already clearly 
evident from the ID FES reported in Figure [9j In the 
case of very small supercells (4x4x6 and 4x4x8 
fee cells, containing 384 and 512 atoms respectively) fi- 
nite size effects arc quite severe and one can hardly see 
a plateau region. The free energy at these supercell sizes 
changes quite rapidly over the whole CV space, probably 
due to the strong interactions between the two interfaces 
formed during the solid-liquid transition, which are quite 
close in such short cells. As we increase the length of 
the supercell at constant xy dimensions, from 4 x 4 x 12 
onwards, another feature of the free energy is observed: 
one can clearly distinguish a plateau with a linear resid- 
ual slope. This can be explained by the reduction in the 
liquid entropy by the constraint of finite xy dimensions 
of the supercell, which slightly raises the melting tem- 
perature. Since the solid is marginally more stable, once 
the interface is formed there will be a linear increase in 
free energy as the fraction of liquid phase grows. Indeed, 
simulations for supercells with larger xy dimensions yield 
a flatter plateau (see Figure [9fb)). The small increase in 
melting temperature is expected to manifest itself when 
the width of the supercell is less than some correlation 
length 2 L corr . L corr can be estimated by looking at 
the distance at which the pair correlation function (see 
Figure [6]) approaches 1, which in our case is 5er, corre- 
sponding to « 3 cell parameters, suggesting we need at 
least 6x6 unit cells in xy. Indeed the effect is seen to 
have vanished by 9 x 9 unit cells (Figure [9Jd). 

With these concerns about finite-size effects in mind, 
we can discuss a reasonable protocol to compute G xs . For 
the 2D metadynamics, the region with 0.35 < sa,sb < 
0.55 is sufficiently flat, and we estimate the free-energy 
of the two-phases configuration to be G s \i — G(sa — 
sb = 0.45). Let us now consider the ID case, where 
region B is restrained to remain solid. In this case, due to 
the finite slope, the estimate of 7^ will depend on which 
point on the plateau we take to be corresponding to the 




Figure 9: (colour online) (a) The plateau region correspond- 
ing to the presence of the solid/liquid interface, with differ- 
ent relative amounts of the two phases, is drawn for different 
supercell sizes in the z direction. The curves are shifted for 
display purpose, and they are only meant to demonstrate how 
the fiat portion of G(s) becomes more extended as larger su- 
percells are considered. Finite-size effects are present also 
in the region for s < 0.1. In inset (b) we compare results 
for simulations with different in-plane sizes, showing a fur- 
ther finite-size effect which depends on the change in T m and 
causes a residual slope of G(sa)/A even after the complete 
formation of an interface. 



value of Gm. Hence, we get our best estimate of G s \i as 
G(sa = 0.2) (the point equidistant from the limits of the 
plateau region) and estimate roughly the systematic error 
due to the finite slope as G{sa = 0.3) — G(sa = 0.1). 

The results of calculations with different supercell 
sizes are reported in Table [TTJ where we report our 
best estimate of j s i and of the statistical and finite-size 
errors(A7|* at and Aj^f respectively). A7|* at is com- 
puted as the root mean square deviation between 10 in- 
dependent simulations of 10 7 timesteps each. The results 
of a run using two CVs is also reported for comparison. 



C. Comparison with other methods 

With the aid of Table |n] we can now discuss the rela- 
tive merits of our technique. First of all it can be seen 
that our calculated value for the (100) surface is very 
close to the ones calculated by CFM(0.369 ± 0.008) and 



10 



# atoms (cell) lsl (AjT'Alii 



51 (2D) 

52 (ID) 

53 (ID) 

54 (ID) 

55 (ID) 

56 (ID) 

57 (ID) 

58 (ID) 

59 (ID) 



2352 (7 x 7 x 12) 
384 (4 x 4 x 6) 
512 (4x4x8) 
768 (4 x 4 x 12) 
1024 (4 x 4 x 16) 
1280 (4 x 4 x 20) 
1200 (5 x 5 x 12) 
2352 (7 x 7 x 12) 
6480 (9 x 9 x 20) 



0.37(0.01) 
0.39(0.0008,0.02 ) 
0.390(0.001,0.008) 
0.390(0.003,0.005 ) 
0.390(0.002,0.006) 
0.390(0.003,0.002) 
0.386(0.003,0.003 ) 
0.369(0.002,0.0006) 
0.360(0.003,0.0004) 



Table II: Value of 7 s j calculated for different supercell sizes 
with both ID and 2D metadynamics. The error is reported 
as (A7f* at ,A7f i 3 ), where A7| ; tat and A7f, s are the statisti- 
cal and systematic error respectively, as defined in the text. 
For all sets of parameters, ten independent runs have been 
performed, each 10 7 steps long. 



CM ((0.371 ±0.003) in RefPand (0.34 ±0.02) in ReP). 
Although we cannot make such a direct comparison with 
CNT (because only an averaged value for 7^, 7™ ff for 
all possible orientations is given) we point out that their 
value of 7™ 9 = 0.302 ± 0.002 is much lower than ours. 
The anisotropy in 7^ accounts for part of the difference, 
as the (110) and (111) surfaces have a lower 7 S P, but 
we suggest that the anisotropy is too small to account 
for all of it. The fact that the value calculated by CNT 
is much lower than both ours and that of the CFM and 
CM may also be due to the curvature and temperature 
dependence of 7^ ; CNT is the only method dealing with 
curved interfaces at temperature below the equilibrium 
T TO , as noted irP. We also point out here that we do not 
neglect the pV term as done in the first version of the CM 
approach^, but we still recover a free energy higher than 
that calculated by CNT. This should rule out the pos- 
sibility, as supposed in Ref.§J, that relaxations in volume 
during the formation of the interface could be another 
explanation for the discrepancies in 7^. 

In part, the existence of a small discrepancy between 
results of CM, CFM and the present work, which rely on 
similar thermodynamic assumptions, can be explained in 
terms of differences in the technical details of the calcu- 
lations. For instance, in some of the CM calculations 
temperature-control has been implemented by a non- 
standard velocity-rescaling method, which might affect 
the accuracy of sampling of the canonical ensemble. In 
the present work we have tried to highlight all the pos- 
sible sources of statistical and systematic error, to facili- 
tate further comparison. In any case, the discrepancy be- 
tween different numerical approaches is negligible when 
compared to the errors affecting experiments, which can 
give results differing by as much as 300% (see e.g. ReP^). 
Hence, any of the aforementioned techniques can be ex- 
tremely valuable in assisting the interpretation of exper- 
imental data and the development of new materials. 

The small system size required for our simulations will 



be a particular advantage, since system size is by far 
the biggest limitation in applying more sophisticated po- 
tentials. We obtain reliable results with system sizes as 
small as about 1000 atoms, more than two order of mag- 
nitude smaller than required by both CFMs and CNT. 
CMs require a few thousand atoms, so the advantage 
is less impressive. However, we remark that the lower 
bound attainable by CM is most likely set by the need to 
mitigate hysteresis effects, while with our metadynamics 
approach this is not an issue, and the limiting factor here 
will be the kind of interactions between interfaces that 
are inevitable in all total energy calculations based on 
periodic boundary conditions. 

In view of the large experimental inaccuracies in the 
measure of 7^; (e.g. Ref!^) even simple empirical poten- 
tials would already lead to a quantitative improvement 
of the knowledge of 7 S /. We are applying our method to 
some of these cases and the results will be published in 
a future paper. Nevertheless, it may be that the accu- 
racy or information given by a self-consistent electronic 
structure method is desired for the interfacial free energy 
calculation, in which case our approach would still be a 
promising candidate. With high performance comput- 
ers, simulating a few hundred atoms for several hundred 
picoseconds is within the reach of present, widely used 
molecular dynamics methods employing so-called ab ini- 
tio (electronic density functional) techniques for the cal- 
culation of interatomic forces. This would probably re- 
sult in better predictive power and smaller overall errors 
despite the possibility of mild finite-size effects. 

With a view to performing calculations with more so- 
phisticated potentials, metadynamics offers a further ad- 
vantage over the other techniques. One could implement 
a process of iterative refining, whereby one performs a se- 
quential set of calculations with potentials of increasing 
sophistication and computational cost, in order to reduce 
the burden of levelling the FES. In fact, the major fea- 
tures of the FES can be captured by the use of very sim- 
ple potentials reproducing the nearest neighbour bonding 
in the real material. This first level FES, G^°^(s), could 
then be used as the initial bias for a second metadynamics 
run, to be performed with a more accurate (and expen- 
sive) potential. At this stage, one will have the much 
easier task of correcting the discrepancy between G {0 \s) 
and the FES of the new potential, G^'(s). This scheme 
could be repeated with increasingly accurate potentials. 



V. CONCLUSIONS 

In the present paper we have presented a novel ap- 
proach to the calculation of the solid-liquid interface free 
energy 7^;, and discussed its application to the calcu- 
lation of 7 S ; for the (100) surface of a Lennard- Jones 
solid in contact with its liquid. Our method is based 
on the definition of a new order parameter, which is de- 
signed to identify /cc-ordering of atoms in the orienta- 
tion of choice, compatible with the periodic boundary 
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conditions, and uses metadynamics simulations to es- 
timate free-energy differences between the bulk phases 
and configurations where macroscopically flat interfaces 
are present. We obtain results for j s i in good agree- 
ment with previously proposed methods. Moreover, our 
technique offers several advantages compared to previ- 
ous approaches discussed in the literature. It requires 
fewer atoms than those methods based on macroscale 
approximations, such as measuring capillary fluctuations 
or the critical nucleation radius, while being less affected 
by hysteresis than cleavage methods, since the interface 
is created and destroyed several times during each simu- 
lation as equiprobablc sampling of the free energy surface 
is approached. We discuss at length the different sources 
of error, and how they can be controlled. In particular, 
we show that our approach is effective even for super- 
cells containing fewer than 1000 atoms, with finite-size 
errors whose importance can be gauged easily. For this 
reason, we speculate that it would be possible to perform 
an ab initio calculation of j s i, at the level of electronic 
density functional theory. To this end, we suggest that 
an iterative refinement scheme, which starts with a bi- 



ased free-energy surface computed from a semi-empirical 
potential, could be a helpful starting point for obtaining 
converged results within reasonable computational time. 
We plan to attempt these calculations in the near future, 
after having further validated our method by comparison 
with experiments in the case of a simple metal for which 
we expect empirical potentials to be adequate. 
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